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I. INTRODUCTION 

Finite amplitude waves in a plasma have been studied intensively for decades in regard to a broad range of physical 
problems related to astrophysics, magnetic and inertial confinement thermonuclear fusion and in nonlinear wave 
theory In particular, nonlinear plasma waves arc of crucial importance for wakefield acceleration in plasma 
configurations where the wakewave is generated either by laser pulses 0, Q or by bunches of relativistic electrons 
[|| , for high-harmonic generation Q and for many other aspects of laser- plasma physics 0, 0] • In order to support a 
strong electric field the Langmuir wave must be highly nonlinear. In a stationary wave the limit on the field amplitude 
is imposed by the wave breaking condition |8[ , while in a nonstationary wave in the regime beyond the wavebreaking 
point the electric field can be even higher Q. 

Nonlinear wave breaking exhibits one of the fundamental phenomena in the mechanics of continuous media. When 
the wave amplitude approaches and/or exceeds the breaking limit the wave form becomes singular as its profile 
steepens, finally leading to the formation of a multi-stream motion. Even in the simplest case of one-dimensional 
electrostatic Langmuir waves in collisionless plasmas this process still attracts great interest due to its importance 
both for the wave amplitude limitation [H, [Tol - [T2| and for its practical relevance to the electron injection into the 
wakefield acceleration phase 0, [l| . In the application to the laser wake field acceleration attention is paid mainly to 
the determination of the upper limit for the electric field d, [IH, [I3 - [l8j . 

Thermal effects in a warm plasma can reduce the maximum wave amplitude 0, Il2l . [I3 - [l8j and modify the character 
of the singularity fl8j . A finite plasma temperature limits the electron density in the breaking wave but in the general 
case does not necessarily lead to smooth density distributions. Since the results obtained by B. Riemann in the 19 th 
century on the wave breaking of nonlinear sound waves (see Ref. [20j . and [H|), it has been known that thermal 
effects do not prevent the "gradient catastrophe" . In this case the singularity in the breaking wave corresponds to 
a shock-like wave profile. Other remarkable singularities are known for nonlinear waves on a water surface which at 
the breaking points become of the type of "Stokes's traveling crested extreme wave" with the interior crest angle of 
2n/3 [22|,[23|]. We also note here the exact solutions, known as "peakons", of nonlinear partial differential equations 
describing the waves on shallow water that have the form of a soliton with a discontinuous first derivative [24| . 

In the present paper we analyze the structure of the Langmuir wave breaking and show that crested Langmuir 
waves in thermal plasmas have a profile with a discontinuous first derivative. 

II. THE WATER-BAG MODEL FOR A RELATIVISTIC LANGMUIR WAVE IN A THERMAL PLASMA 

A. Electron distribution function formed as a result of a gas multiphoton ionization 

In the case of a plasma irradiated by a high intensity laser pulse the temperature is determined by the laser light 
parameters for the time interval before the main pulse comes. During the interaction of a femtosecond, tcrawatt 
laser pulse with gas targets a plasma is created via photoionization [25( by the prepulse or by the ASE (Amplified 
Spontaneous Emission) pedestal. In such a collisionless plasma the electron energy is of the order of the quiver energy 
in the ionizing laser field, i.e. typically in the range below keV. It is thus substantially lower than the electron energy 
in the main laser pulse, which is typically in the MeV range, that excites the wake plasma wave (e.g., see Fig. [TJ where 
a typical electron distribution function formed as a result of optical field ionization of the gas target by an ultrashort 
laser pulse is shown, (26j). Being limited by the quiver energy, the electron distribution function is not Maxwcllian 
and can be adequately described by a simple water-bag model, which is otherwise considered to be too artificial and 
restrictive. We note that the water-bag electron distribution function has been used in Refs. 

B. Basic equations 

Following Ref. , we consider the electron phase space {x,p) shown in Fig. [2 which corresponds to the support 
of the electron distribution function / e (p, x,t). 
The electron distribution function is constant 

f e (p, x, t) — constant (1) 

within the region with borders marked by p+(x,t) and P—(x, t), while f e (p,x,t) = outside this region. Here the 
constant is proportional to the ratio of the electron density and the momentum width. The electron distribution 
function can also be expressed via the unit step Heaviside functions: 



f e (p, x, t) = constant x 9(p — p~(x, t)) 9(p+(x, t) — p), 



(2) 
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FIG. 1. Electron distribution function f e {Pz) (arb. units) formed as a result of optical field ionization of the gas target by 
a 40 fs laser pulse with dimensionless amplitude eE/m e uJoc = 0.1, which is the tunneling regime of yx <C 1 [25|]. The f e 
dependence is on the normalised electron momentum p z /m e c along the laser polarization direction. 



where 0(x) = for x < 1 and 6(x) = 1 for x > 1. 

The evolution of the distribution function is described by the Vlasov-Poisson system of equations 

dtfe + vd x f e -Ed p f e = 0, (3) 



d x E = 1 - n e , (4) 

where all the variables are written in a dimensionless form normalised in a standard way in which the time and space 
units are uj~^ and cw"/, the momentum and the velocity are normalised on m e c and c, the unit for the electric field, 
E(x, t), is m e oj pe c/e, with w pe = (47rnoe 2 /m e ) 1 / 2 being the Langmuir frequency, e and m e are the electron charge and 
mass, and no is the density of ions which are assumed to be at rest. The electron velocity is equal to v = p/ (1 +P 2 ) 1 / 2 
, and n e (x,t) is the electron density normalised on no. Global charge neutrality is assumed. Eq. ([3]) describes the 
incompressible motion of the distribution f e in phase space. 




FIG. 2. Electron phase space in a single water-bag model. The electron distribution function is constant f e (p,x,t) = g e within 
the region with borders marked by p+(x,t) and p-(x,t), while f e (p,x,t) = outside this region. 



Calculating the first momentum of the distribution function we find that the electron density is related to the 
bounding curves p+(x, t) and p_ (x, t) as 

n e (x,t)= J f e (p,x,t)dp = g e \p+(x,t) - p-(x,t)], (5) 

— oc 

where g e is a numerical constant (see Eq.([T])) that gives the ratio between the dimensionless electron density n e (x,t) 
and the dimensionless momentum width Ap(x,t) = p + (x,t) — p-(x,t), and is determined by the value of this ratio 
at t = 0. 
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From Eqs. (|3]) and ((4]) it follows that the functions p+(x, t), p-(x,t) and E(x,t) evolve according to (sec also 
Ref.[l(l) 

d tP+ + 1± d x p+ = -E, (6) 



dtp- + j£ d xP - = -E, (7) 



d x E=l-g e (p + -p_). (8) 

C. Dispersion equation for the wave frequency and wave number 

A large energy spread of the electron distribution function leads to a change of the Langmuir wave frequency due to 
its dependence on the plasma temperature. Linearization of Eqs. ([6] [8]) around the equilibrium solution p + = p+,o, 
P- = p-,o, E = with p-,o = — P+,o = Apo/2 gives the dispersion equation for the frequency, u>, and the wave-number 
k in the case of the small amplitude Langmuir wave. In dimensional form it can be written as 

^ 2 = 87rn e 2 c k 2 c 2 Ap 2 

y/Amlc 2 + Ap 2 4m2 c 2 + Ap 2 ' 

The corresponding kinetic dispersion relation for a relativistic Maxwellian distribution function (Jiittner-Synge dis- 
tribution) and its derivation in terms of relativistic fluid-like equations are given in Ref. [13] and references quoted 
therein. In particular waves with phase velocities larger than the speed of light are considered in Ref. [13] since in the 
case of the Jiittner-Synge distribution, and in general of a distribution that is not piece-wise constant in momentum 
space, waves with phase velocity smaller that the speed of light are heavily damped in the relativistic regime [28| . 
It is shown that for these "superluminal" waves the relativistic electron population obeys an effective isothermal 
equation of state as soon as the normalized thermal momentum p t h becomes larger than an appropriately redefined 
phase momentum p p h = (jj/{lu 2 — k 2 c 2 y/ 2 . We note that although the superluminal regime could also be considered 
within the water-bag formalism by redefining 7 P /, below Eq. (|22]) . such a regime is not of interest for the investigation 
of wavebreaking. In fact, both Landau damping and wavebreaking are related to particles that have an unperturbed 
velocity (in the case of the linear Landau damping), or are accelerated by the wave electric field to a velocity that 
matches the wave phase velocity and this cannot occur for superluminal waves. Conversely, this indicates that ther- 
mal effects tend to favour Landau damping in the case of a particle distribution that is not piece-wise constant and 
wavebreaking in the case of a water-bag distribution. 

As can be seen from Eq.(|9]) a finite temperature modifies the Langmuir frequency and makes it depend on the 
wavenumber, k. In terms of the variable X = x — v p ht the wave is characterised by the wavenumber k w , which is 
given in dimensionless form by 



kw — , / O o 9 2 ' (10) 

where 7+^ = y/l + Pq- The wave number k w tends to infinity for p+,o/7+.o — > /?ph ■ The frequency dependence on 
the electron temperature leads to a shortening of the wakewave wavelength. It results in a lower electric field in the 
wakewave in comparison to the case with a relatively small thermal spread. 

D. Long wavelength limit 

It is easy to obtain from Eqs. ([6]-|8]) that spatially homogeneous nonlinear oscillations of electrons in relativistic 
thermal plasmas are described by the system of ordinary differential equations 
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dE _ ^ l +P\- V 1 + (P+- A Po) 2 
dt Apo 

where Apo = p+fl — P—fl. The Hamilton function corresponding to these equations is 



(12) 



H(E, P+ ) = ^-+TL( P+ ) 



with the potential function 



where the function U (z) is given by 



n( P 4 



u(p+) - u( P+ - A Po ) 

2A Po 



U(z) = z\f\ + z 2 + In (z + \/l + z 2 ^j . 



(13) 



(14) 



(15) 



Isocontours of the Hamiltonian function ([13)) in the plane E,p + are shown in Fig. for Apo = 3. The potential 
function II (p + ) is plotted in Fig. [3]}. Nonlinear oscillations are shown in Fig. [3J;, where the time dependence of the 
electron momentum, p+(t) and electric field, E{t) are plotted for p+,o = 12.5 and Apo = 7.5. The momentum p + 
oscillates between the value p + $ and — (p+.o — Ap )- 
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FIG. 3. a) Isocontours of the Hamiltonian function (|13[) in the plane (E,p+) for Apo = 3; b) Potential function Il(p+; Apo); c) 
Time dependence of the electron momentum, p+(i) and electric field, E(t), for p +i o = 12.5 and Apo = 7.5. 



In the small amplitude limit the oscillation frequency is given by expression Q for k = 0, i.e. 



(87rnoe 2 c) 1 



(4m2 c 2+Apg) 1/4 ' 

If the oscillation amplitude, p m , is large compared to the half- width of the distribution function, Apo/2, and p m 1 
then the oscillation period is given by T = (4p m /7m e 2 ) 1/2 (see Ref. @). 



E. Langmuir waves travelling with constant velocity 

We consider waves propagating along the x axis with constant phase velocity v p h, where all functions depend on 
the independent variable 



In this case Eqs. (H3 -[5J take the form 



X = x - (3 p ht. 



h'+ = -E, 



(16) 
(17) 



ti = -E, 



(18) 



and 



(19) 



where we introduced the dependent variables h + and /i_ defined by 



h±(p±) = y 1 +p± - /3pkP±, 
a "prime" denotes differentiation with respect to X and /3 p h = v p h/c. We use 

Apo =P+,o -P-,o, 

with p + ,o = p + (Xq) and p-.o = P-(Xq) taken at X = Xq, where E 1 = 0. Inverting Eq. ([20)) we obtain 



p±(h±) = 7ph/3ph/i± - Jlph h ± ~ 7ph> 



With 7ph = 1/y/l-Pl 

Here and below we assume "subluminal" propagation velocity, i.e. /3 p h < 1. 



F. Hamiltonian form of the equations describing a travelling Langmuir wave 
Multiplying Eq. (|19|) by E and using Eqs. (jTTJ) and (fT8|) and integrating it over X, we obtain the integ 

^ + k -±(l- + k -(l + ^f hh -) + " = constant, 

where 

= z^z 2 - 1 - In (z + Vz 2 - 1) . 
This function vanishes at z = 1. In the limit z — > 1 + its behaviour is described as 

W{z) = \V2{z - 1) 3/2 + 0{z - 1) 5/2 . 

For 2 — > oo we have 

W{z) = z 2 -\-\n2z+^ + o(± 

According to Eqs. (fT7|) and (fT8|) the variables and h + , are not independent and are related by 

h- = h+ — A/io, 

where the constant Ah Q is determined by the values of p + and p_ at X = X a : 



Ah = h +fi - h-fi= ^Jl+p 2 +Q - ^Jl+p 2 _ - /3 p h(p+,o - P-,o)- 
As a result, Eqs. (|17[ [PS]) and (|T9| can be rewritten in the form 

h+ = -E, 



E i = l_ P+( h +) - P-( h + - Ah o) 
Ap 



This is a Hamiltonian system with Hamilton function 



n(E,h + ) = ^- + u(h + ), 
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where h + and — E are canonical variables, and 

L_ 7g h ^\ + ^^-^(^-AM). (32) 
V Ap y 2Ap 

For a symmetrical distribution where p+,o = —p—,o Eq. ([28| takes the simpler form Aho = — /3 p hApo and the 
potential H(h+) reduces to 

n (h ^ ji . , ^(7 P hfe+) - ^(7 P h(fe+ + IphApo)) , . 

ll S ym(, ft +) = 7ph / M ■ (33) 

The potential n(/i + ) is plotted in Fig. 0^, as a function of h+, /3 p h for four values of Apo assuming that p + .o = —p-.o, 
i.e. Ah = — /3phApo- Isocontours of the Hamiltonian function in the plane E,h + for /3 p h = 0.8 and Ap = 0.1 are 
shown in Fig. 0b. 




III. THE WAVEBREAKING LIMITS 
A. Crested Langmuir wave 
The system of Eqs. (|17H19I) has a singular solution when h+ — » 7^, i.e. 

P+ -> P+m = I ^ Ph (34) 

which corresponds to the wavebreak in thermal plasmas when the electron velocity calculated for the momentum on 
the upper bound curve, p+(x, t), becomes equal to the wave phase velocity. In this limit dh + /dp + and the upper 
bound curve is no longer a single valued function of X . 

The electron momentum on the lower bound curve at wavebreak is 



P-M = P+m - /3 P h7phA/Jo - V /7ph A ^ - 2 7 ^ h A/i . (35) 
Eqs. ([5l [34] [35]) give for the electron density at the wavebreaking point 

n e ,br = t ■ (36) 

Ap 

For a symmetric distribution function such that p+,o = — P-.o from the electron density dependence on h, 

7 P h 



T^(ft + /3p h Apo) a -l-^fc 2 -l 
Apo" 



n e (h) = LV — : L - 7ph£ph, (37) 
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it follows that at the wavebreaking point, h — > 7 ph 1 ! the density tends to 



n e>br =7 P V ph Ul+^-^-fa . (38) 



In the nonrelativistic limit, when /3 p h <C 1, Apo *C l//3ph, and 7 P h ~ 1, the density is 



In the ultrarelativistic limit, when /3 p h ~ 1, and 7 P h > 1 we have 



n £ M « W ^rr^ - 7ph^ph- (39) 



provided Apo <C 2//3 p h7 P h (see also [H[) while for Apo = 2/3 p h7 P h we have n e ,br = lj because in this limit a wave with 
arbitrarily small amplitude breaks, as seen from Eq. (fTU|). In the above considered cases the electron density written 
in dimensional units is 



»e,br ~ n . I " 2eVph (41) 
P+,0 



for a nonrelativistic plasma wave and 



n e ,br«noy (42) 

V P+,o 

in the limit 7 P h 3> 1. 

In order to find the density behaviour in the neighbourhood of the breaking point, we expand the electron momen- 
tum, on the upper bound curve, in the vicinity of its maximum, SX = X — X^ T — > 0. Here X^ is the location of 
the breaking point. Locally, the momentum is represented by 

p+=p+M-Sp++0(Spl) (43) 

with p+,br given by Eq. (fM)) . Keeping the main terms of the expansion over 8p+ of Eq. (j3~?| we obtain for the electron 
density 

n e = ra e ,br - (44) 
Ap 

where we used the expression h tn 7^ + <5p+/27 ph . From Eqs. (|29 [ [30 ]) for the dependence of Sp+ on X we have 

(<$p* )" = 2n e>b r7ph- (45) 

Integrating this expression we find 



Sp+=± x /n eM r^SX, (46) 



where we assumed that 5p' + at SX vanishes, i.e. the electric field at the breaking point is equal to zero. Since by 
assumption 6p+ must be non-negative, we must chose the " — " sign in the interval SX < and the " + " sign for 
SX > 0. As a result we can write for the momentum p+ in the vicinity of the wavebreaking point 



P+ ~ P+,b 



r" ^e,W ph \SX\. (47) 



From the expression ([5|) for the density, and recalling that at wavebreak <5p_ cx (Sp + ) 2 , we find that in the vicinity of 
the breaking point the electron density can be written as (see also (T^ |) 



"1 e ,br7p h 

He ~ n e ,br - |^|. (48) 

Apo 
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FIG. 5. Structure of the nonlinear wake wave: a) electron phase space, b) the electron density (in the inset the density 
distribution is shown in the vicinity of the maximum), c) electric field as functions of the coordinate X. The normalised wave 
phase velocity is /3 p h = 0.999; the plasma thermal momentum width at X = 15 is Apo = P+,o — P-fi = 5; the maximum of the 
electric field is E m = 0.4190625. 



This type of wave breaking in the general case corresponds to the " peakon" structures known in water waves [22h24| . 
It can also be called "A- type" breaking. 

The structure of the nonlinear wake wave of the electron density, the electron phase space and the electric field 
is shown in Fig. [5j as obtained by numerical integration of Eqs. ([6]-|8|). In Fig. [5] we show the high temperature 
case with the initial distribution function width p+,o comparable with the value of electron momentum on the upper 
bound curve at the wavebrcaking point, p^ T . From Fig. [SJd we see that in the vicinity of the density maximum (see 
inset to Fig. the density dependence on X corresponds to Eq. (|48l) . 



B. Maximum electric field in stationary wave 



As is seen from the trajectory pattern in the E, /i + plane presented in Fig. 0}d, the electric field maximum is 
reached at the point h + = He (see also Fig. [5]) where the derivative of the electric field with respect to h+ vanishes, 
dE/dh+\ h _ hE = 0. This condition results in the equation for He'- 

Ap =p + (h E ) - P+ (h E + f3 ph Ap ), (49) 

where the function p+(h E ) is given by Eq. (|22l) . Here we assume the symmetric distribution with p+,o = — P-,o = 
Ap /2. The solution of Eq. (gSJ) is 

h E = fi^M-feg*. (50) 
The last bound trajectory in the E 1 h + plane is determined by the equation 



^ + n(^ + ) = n( 7p - 1 ) 



(51) 
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FIG. 6. a) Isocontours of the Hamiltonian function (Hi = 1.9125 and 2.25) in the plane (E, h + ) for /3 p h = 0.6 and Apo = 0.05. 
The electric field maxima -E max and -E maXj AB are reached at h+ — He- b) The curves 1,2,3 indicate the potential 11(h) 
values corresponding to (1) a periodic wave with amplitude below the wave breaking limit, (2) to a wave at the wave breaking 
threshold, and (3) to a wave with amplitude above the wave breaking limit. Vertical dashed line marks the wave breaking 
boundary. 



with the potential given by Eq. (|33j). Substituting /i+ = He we find the electric field maximum 



n( 7ph 1 ) - u(h E ) 



(52) 



In the limit of cold plasma, when Apa — > 0, i.e. p± — >■ p with p±_o — > the Hamilton function pip reduces to 

H(E, h) = ^- + 7> - fay/h?^ ~ 7* h > (53) 



where h = y/1 + p 2 — f3 p hp. It can be rewritten in the form of the energy integral, E 2 /2 + 7 = constant. 

The potential H(h) = 7^/1 — f3 p h\/ ^ 2 7ph ~~ 7ph * s plotted m Fig- as a function of h and /3 p h- Isocontours of the 



Hamiltonian function in the plane E, h for ,9 p h = 0.5 are shown in Fig. [7t>. 
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FIG. 7. a) Potential function H(h) vs h,; b) Isocontours of the Hamiltonian function in the plane (E,h) for f3 p h = 0.5. 



It is easy to see that the electric field equals zero at the maximum of the electron quiver energy. This condition 
yields the result obtained in Ref. for the maximum value of the electric field: 



E AP = y2( 7ph - 1). (54) 
For small but finite electron temperature, Ap <C l/(/3 p h7ph)j we obtain (see also Ref. fTij ) 

g max « J2(~f ph - 1) - \ ^^ 3/2 ^Ap~ . (55) 
v 3 V(7 P h - 1) 

At ApQ — > 2/3 p h7 P h the electric field vanishes because, as mentioned before, in this limit the wave with arbitrarily 
small amplitude breaks. 
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In Fig. [8] we show the maximum electric field, E max , in the breaking wake wave. The dependence of this field, 
normalized on -Eapj on the wave phase velocity /3 p h and on the width of the electron distribution function Ap is 
presented in Fig. [8] a. Figures [8] b and c show dependences of E max and Eap on 7 p h for small and large Apo- We see 
that in the limit 7 P h 3> 1 the difference between E max and Eap increases according to Eq. (|55|) . Figures [5] a and c 
clearly illustrate the above mentioned fact that in the limit Apo — > 2/3 p h7 P h the value of E max vanishes. 




FIG. 8. Maximum electric field in the breaking wake wave: a)_E max normalized on Eap, depending on the wave phase velocity 
/3 p h and the width of the electron distribution function, Apo; The curve in the /3 p h - Apo plane where E max = is given by 
Ap = 2^p h 7 ph . b) £J max and Eap v.s. 7 ph for Ap = 0.125. c) _B max and Eap v.s. j ph for Ap = 7.5. 

Independently of whether the plasma temperature is finite or vanishes, from Eqs. (|3 1 [) and (foU]) we obtain that 
the second derivative of the potential II with respect to h+ (h) becomes singular at h + = 7^ (h = 7^ ), which 
corresponds to the vertical (dashed) singular line in Figs. [3Jd and[B}D. In this limit the Hamiltonian in Eq. ([53]) takes 
the value 

E 2 

— +7 P h, while 7 (/i = 7^) = 7 ph . (56) 



C. Cold wavebreaking limit 



In order to compare the properties of the singularities formed in thermal and cold plasmas we reproduce here the 
dependence of the electron momentum and density on the coordinates in the cold wavebreaking case (for details see 
Rcf. |19j). In the cold plasma with p +0 and p_ 0> which implies p+(X) = P-{X) = p, equations (fTTl — fT^|) 
can be reduced to 

P 



1 + p2 -^> ' ^VTTF-p - (57) 

The solution of this equation can be expressed in terms of elliptic integrals. In order to analyze these solutions in the 
vicinity of the singularity we note that its right-hand side becomes singular when the denominator, /? p h \J 1 + P 2 — P, 
tends to zero, i.e., when the electron velocity v becomes equal to the phase velocity of the wake wave. In the wake 
wave, the singularity is reached at the maximum value of the electron momentum, p m = p^ r . We assume that the 
singularity is located at the coordinate X = X^ r . We consider the wave structure in the vicinity of the singularity, 
and find that here the electron momentum depends on SX = X — X^ as 

P = P^-^lUi\6X\/2^ 2 ) 2 /\ (58) 

The electron density tends to infinity as 

n » 7io7ph(2 1/2 /W3|<5X|)- 2 / 3 . (59) 

The 2/3 power behaviour can be recognized in Fig. [9J which presents the wakewave generated in a relatively low 
temperature plasma with p+ ; o ^ Pbr- However, from Fig. ^jp we see that in the very vicinity of the density maximum 
(see inset to Fig. [5b), the dependence of the electron density and momentum on X still corresponds to Eq. 
showing at the wave crest the density profile which can be approximated by the "peakon" dependence. The electron 
distribution width, p+,br — P-,t>n near the maximum is characterized by the value 



Ap = p+,br - P-,br ~ y 2/? p h7 ph Ap , (60) 

where we assumed 7 p bApo <C 1. 
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FIG. 9. Structure of nonlinear wake wave: a) electron phase space, b) the electron density (in the inset the density distribution 
is shown in the vicinity of the maximum), c) electric field as functions of the coordinate X. The normalised wave phase velocity 
is /3 p h = 0.995; the plasma thermal momentum width at X — 15 is Apo = P+,o —P-.o = 0.25; the maximum of the electric field 
is Em = 1.9876. 



IV. HYDRODYNAMIC APPROACH 

A. Waterbag distribution 

The system of Eqs. ([6] -[8]) can be written as a system of hydrodynamic-type equations 

d t N + oU = 0, 

ftP + d x G = -E, 

d,E + ¥d x E = V 



for the electron density 



average momentum 



and electric field E. Here 



N(x,t) = J f e (p,x,t)dp = -^—(p + - p-), 



(P+ +P-), 



J(a,*)= / 7= „ fe(p,x,t)dp = — i— (7 + -7-): 



(61) 
(62) 
(63) 

(64) 

(65) 
(66) 
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with 7 ± = Jl+p±, 



and 



Y(x,t) =1 / _^ = ^/ e ( P)a:>t )dp = ^i±^ ) (67) 
NJ Vl+P 2 7+ +7- 



Pfe(P,X,t)dp 

fe(p,X,t)dp 



— v = i( 7+ +7_), (68) 



These functions arc related to each other as 



and 



(69) 



In the case of a wave travelling with constant velocity c/3 p h, the functions N, P, G and E depend on the variable 
X = x — fipht and Eqs. ([6T1- l63"| can be reduced to 

c-'W— S^p- (71) 

N /? ph G-P' (?2) 

where a prime denotes differentiation with respect to X. Equation (|7T|) looks identical to Eq. (|57|) which describes 
the wave break at p/ 7 — > /3 p h and the formation of a singularity in the electron density, n — > 00, according to Eq. 
(fBT))) . However, due to the nonlinear dependence of G on P given by relationships ([70)) and (|72fl the character of the 
singularity changes and becomes of the type described in Sec. MI Al In particular, we can see that the condition for 
the denominator in the r.h.s. of Eq. (|7ip to vanish implies that V = /3 p h- This condition can be rewritten as 

A. = P - ±± ^~ = ^ ( 7 ± ^ ± ) • (73) 



7+ + 7- 7+ V 1 +7-/7 

Assuming that in this limit p_ = p + + Sp with 5p/p+ <C 1, we can easily find that the condition of " hydrodynamic 
type wave break" ([73"]) used in [TJJ is equivalent to 

^ = ^ = ^(l-%), (74) 
7+ +7- 7+ V P+J 

which requires > /3 p h, i.e. the waterbag description in the adopted limit of a stationary nonlinear wave 

propagating with constant velocity is no longer valid. 

B. Nonrelativistic limit 

In the nonrelativistic limit Eqs. ([61]) and (|62|) take the form (see Ref. (Iol |) 

d t n + d x (NV) = 0, (75) 

An 2 

d t V + Yd x Y = -E- -Jp-d x N 2 , (76) 
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which corresponds to a gasdynamics system where the pressure depends on the gas density as 

with Po=const. 

For a wave travelling with constant velocity, /3 p h, we obtain 



(77) 



- /? P h)~ - 



8(V-/3 ph ) 2 



V 



/3 P h-V 



The singular points of this equation correspond to 



h = /3 P h and V 2 , 3 = /3 P h ± 



/3phAp 



(78) 



(79) 



We see that the points Vi and V2 lay beyond the applicability range of the waterbag model, while a wave with 
V max — > V3 is qualitatively described by Fig. [9] with maximum density n max = y / 2 / S p h/Apo and electric field 

-Bmax « /3ph - \//3phApo/2. 

C. Ultrarelativistic limit 

This case corresponds to the limit P> 1. Expanding Eq. (|70]) into series of the small parameter NApo/P <C 1 we 
obtain 



8(1 + P 2 )3/2- 

Using Eq. ([72} for the electron density N we find from Eqs. {71]) and ([SO]) 



Vl + F2-/3 P hP + 



8v ^iP2( / 3 ph yr 



The singular points of this equation, written in terms of the average velocity V = P/Vl + P 2 , are given by 



r i « /3ph and V 2 , 3 « /3 P h ± 



/ /? ph Ap 



(80) 



(81) 



(82) 



with V3 corresponding to the wake wave breaking. This yields for maximum density n max y'2/3 p h7 ph /Apo and for 
the electric field E max ss \/2(7 P h — 1) — (2/3)7 P h-\//3 p hApo m agreement with Eqs. P0|) and (|55j). 



V. COMPUTER SIMULATION OF THE PLASMA WAVE BREAKING IN THERMAL PLASMAS 



During the irradiation of underdense plasma targets by high-power laser pulses, the light within the pulse generates a 
finite amplitude wake wave whose parameters depend, in particular, on the plasma temperature and on the interaction 
geometry. A thorough study of these effects require computer simulations. We performed parametric studies of the 
laser pulse interaction with underdense targets using a two-dimensional (2D3P) particle-in-cell (PIC) code [29| . 

Here the effects of the finite electron temperature have been taken into account in three limiting cases. In the first 
case the initial electron temperature has been assumed to be equal to zero. In the second case thermal effects have been 
modelled by the electron distribution corresponding to the initial waterbag distribution function with a temperature 
equal to lOOcV. In the third case, the initial electron distribution was Maxwellian with the same temperature. In both 
the cases of the waterbag and Maxwellian distributions, the total electron energy is the same, i.e., average energy for 
the waterbag case, 



(mec2 _ ^ / JP_\ = rV/^eW = SArf 

V U }> \2m e / J Ap °P 2 d P 40m e ' 



(83) 



is set to be equal to that in the case of a Maxwellian distribution, < p 2 /2m e >= (3/2)kgT. Here we assumed that 
p/mc <C 1. 
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FIG. 10. Results of 2D-PIC simulations of the nonlinear wake wave generation in thermal plasmas: electron density and phase 
space in a cold plasma (a,d,g); in the plasma described by the waterbag distribution with a temperature of 100 eV (b,e,h); in 
the plasma described by the Maxwellian distribution with a temperature of 100 eV (c,f,i). 



In these simulations, the laser pulse has a normalized amplitude of ciq = eEo/m e ujc = 4.6, a wavelength of A = 
2ttc/lo = 0.8 /xm, focused onto a spot of the size of 13/xm, and duration of 16 fs. The plasma density equals 4x 10 19 cm~ 3 . 
The width of the simulation box is equal to 50 x 65 A 2 . The mesh size is Aa; = A/160 with 30 particles per cell. 

Simulation results for the parameters of interest are shown in Fig. 1101 Here the x coordinate is measured in 1 /im, the 
electron momentum p is normalized on m e c, and the density is normalized on the critical density n cr = m e uj 2 /Aire 2 . 
The figures are plotted for the time when the highest density is reached, which is 350 fs, 310fs, 310fs for zero- 
temperature, waterbag, and Maxwellian distribution, respectively. In the cold plasma case, the electron density 
distribution in the first maximum of the breaking wake wave takes a cusp- like form (see Fig. [10] a. In Figs. [10] b and 
c we see that the finite temperature effects lead to a decrease of the maximum electron density in the breaking wake 
wave, to the broadening of the maximum and to the formation of peakon-like structures for both the waterbag and 
the Maxwellian distributions. Note here more efficient electron injection in the finite temperature plasma compared 
with the cold plasma case. 



VI. ABOVE THE WAVEBREAKING LIMIT 



A. Maximal electric field 



The limiting electric field given by Eq. (|54"|) corresponds to a stationary Langmuir wave for which the electron 
quiver energy is below m e c 2 7 P h- When the Langmuir wave is excited by a short laser pulse its amplitude and its phase 
velocity depend on the plasma density and on the laser pulse intensity Propagating in an underdense plasma, 
an intense laser pulse can accelerate plasma electrons longitudinally up to the energy m e c 2 (l + a 2 /2). In a cold 
plasma the wavebreaking condition corresponds to (1 + a 2 /2) = 7 p h, where 7 p h = (n cr /no)(l + a 2 ) 1 / 4 is the Lorentz 
gamma-factor calculated for the wake wave phase velocity which is equal to the laser pulse group velocity. Here the 
dependence of the electromagnetic wave group velocity on its amplitude is taken into account according to §. This 
yields the wake wave breaking threshold in terms of the driver laser pulse amplitude [3(| 



1/3 



a > 



n 



(84) 
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In general, the laser pulse amplitude in a plasma is different from its value in vacuum due to the laser pulse self- 
focusing and self-channelling [3l[ . The laser pulse amplitude inside the self- focusing channel relates to the laser power 
V and plasma density as |33| 



> 8 t5~~ — 



(85) 



where V c = 2m 2 c 5 /e 2 



i 17GW. 

Using Eqs. (|84|) and (|85|) we obtain the wake wave breaking threshold: 



Zir \ no 



(86) 



For example, from Eqs. 



and 



we find that if a laser pulse of the wavelength A = 0.8/^m, for which 



2 x 



propagates in a plasma with density uq = 2 x 10 cm 



, the wavebreaking threshold is reached 
10 19 W/cm 2 . 



for V > 30 TW and ao = 7.3, i.e. for a laser intensity of the order of 7 x 

A laser pulse with power larger than that given by the r.h.s of Eq. (|86p . causes the wake wave to break in the first 
period, with the electric field well above the limiting value given by Eq. ([5"1[) and with a number of electrons piled 
up in the singularity region much larger than in the stationary case described by Eqs. ([36]) and (|59|) . This fact has 
important consequences for determining the laser wakeficld acceleration scaling [2, HJ • 

A wakewave with an amplitude above the wave break threshold is transient and forms a region with multi-stream 
electron motion. The multi-stream motion region expands in the forward direction at a relative velocity dX/dt « 
c(l — /3 p h) ~ c/2j 2 h . Since in the limit 7 P h 3> 1 this velocity is low, the region with a large electric field (and with 
a large number of electrons) can exist for a substantially long time, which is of the order of the charged particle 
acceleration time, i acc = 2X w j 2 h /c. Here X w is the wake wave wavelength. 

The structure of the wake wave both below and above the wavebreaking limit can be revealed from the phase plane 
pattern presented in Figs. 0] b and [7] b. The stationary (periodic) waves correspond to the bound trajectories in 
the phase plane shown in Figs. 0]b and b. The stationary breaking wave is described by a last closed trajectory 
touching the vertical (dashed) singular line, corresponding to h+ = 7p h j in this figure. 

In the vicinity of the singular line in Fig. |4Jd the Hamiltonian function ([31 [I with the potential in the form given by 



Eq. (p?3"[) can be expanded in series of Sh + = h + — 7 ph as 



H(E,h) 



2 



ni 



~dh 



5h 



(87) 



Here we assume a symmetric electron distribution at X = Xq, where p+,o = ~ P-,o- 
For a finite temperature plasma in the vicinity of the singularity we find 



P h 



2/3 P h 



1/2 



- 1 



5h+, 



where the constant term 7 p h — W(l + /3 p h7 P hApo)/(2Apo) has been dropped. At the wavebreaking threshold, the 
value of the Hamiltonian H = H i = and the electric field E tends to zero at 5h + — > +0 as 



E 



7 P h 



P h 



2/3 ph 
7phAp 



1/2 



1/2 



y/2Sh 



+ • 



(89) 



For Hi > the electric field E\ at wave break, h+ = 7~ h , is given by E\ = \J2H\. 

The quantity of the Hamiltonian Hi is determined by the parameters of the laser pulse driver generating the 
wake wave. In the limit of large laser amplitude, n > 1, assuming that the laser pulse has an optimal duration, 
■Has ~ a/2w po , we can find that Hi fts 4irnom e c 2 a 2 , i.e. the maximal electric field is given by E max w a\/ 8Trnom e c 2 . 

For Hi < 0, the wave breaking condition is not reached, and the electric field vanishes at 



5h+ = 5h + i 



Hi 



7 2 h [^ h + (2/3 P h)/(7 P hApo)] 1 /2- 7 2 h ] 



as 



E = 7 ph 



2 2^ ph 



1/2 



1 1/2 



(2Sh-26h + i) 1/2 . 



(90) 



(91) 
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In the cold plasma limit the Hamiltonian (|53[) expansion in the vicinity of the singularity has a different behaviour: 

H(E, h) = ^-- fth7p( 2 V2a, (92) 

where the constant 7 p h term has been dropped. At the wavebreaking threshold, the value of the Hamiltonian % = 
Hi = and the electric field E tends to zero at Sh +0 as 

E = [/? ph (2 7 ph) 3/2 ] 1/2 (^) 1/4 , (93) 

For Hi > the electric field Ei at wave break, h = j~ h , is given by E\ = \/2Hi. For Hi < 0, the wave breaking 
condition is not reached, and the electric field vanishes at Sh = Shi = E 2 /(2f3 2 h ^ h ) as 

E = [/3 P h(2 7 ph) 3/2 ] 1/2 (V^- VShi) 1 ' 2 . (94) 

In the limit of a relatively low plasma temperature Apo "C /3ph7ph in order to estimate the maximum electric 
field we can use the Hamiltonian in the form given by Eq. (|53[) . In this limit the maximum electric field, £ max = 
\J2 (n(/i m ) — n(ft,E)) with /ie = 1 and h m determined by the maximal electron quiver energy in the wake (see Fig. 
[6]), is given by 



V2 (7m -1). (95) 



The electric field at the wake wave breaking point is equal to y 2 \Jl(h m ) — n(7 ph 1 )j , which in the limit of a relatively 
low plasma temperature yields 

E hl = v /2( 7m - 7ph ). (96) 

For a wake wave with a large enough amplitude, when 7m 3> 7p h, both the maximum electric field and the electric 
field at the breaking point can be substantially larger than the electric field in the stationary wake wave given by Eq. 

As we see in Figs. [4] and [7] in the regime under the consideration the injected electrons appear in the region 
h > hbr = 7p h with a large accelerating electric field. 

At wave break electrons are injected into the region h > hbr = 7 P h 1 where there is a large accelerating electric 
field, as seen in Figs. @] and [7] Note that since the electric field at the breaking point does not vanish, the type of 
the singularity that is formed in the electron momentum and density distributions changes. In a finite temperature 
plasma the electron density in the vicinity of the singular point is determined by Eqs. (|44[) and (|45|). From Eq. (|45|) 
we find 

S P+ = -\/7>e,br<5X 2 + 2 7 3 £ br <Df (97) 

where E^ Y is given by Eq. (f9"6")l and it is assumed that 5X > 0. Inserting Eq. (|97p into Eq. (l4"4l we find that for 
-EW 7^ the electron density near the singularity behaves for SX —> +0 as 

n e = n e ,br - ~t s/2 1 l l E hr 5X. (98) 

In the limit of cold plasma, Ap — > 0, Eq. (|57|) yields (see also [l9j ) 

(6p 2 r = -^f^. (99) 
op 

Multiplying the left- and right-hand sides of this equation by (Sp 2 )' and integrating over X, we obtain 



6 P Sp' = yjTr^Ete - 2p phl * h Sp. (100) 

For Eb r ^= the main term in the expansion of the solution of Eq. (|100[) for 5X — > +0 is 

6p=-(8^ h E hl ) 1 ^V6X. (101) 
Using this relationship we find that in the vicinity of the singularity the density depends on SX as 

- - ^ (102) 



(8^ ph E hl )^VSX 

If instead E^ = 0, the electron momentum and density are given by Eqs. (|4"T1 14"5)) for Ap n ^ and by Eqs. 
for Apo = 0, respectively. 
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B. Results of simulations with the 1-D Vlasov code 



The Vlasov- Poisson system is solved for the electron distribution function, f e (x,v,t), with the numerical scheme 
described in Ref. [32], limiting our study to the 1D-1V case. The equations are normalized by using the following 
characteristic quantities: the charge e and the electron mass m e . The electron density is normalized on the density 
of ions no, which are assumed to be at rest. Time and space coordinate are normalized on the inverse Langmuir 
frequency u>~^ and on the Debye length \d = \/T e /ATm e e 2 , respectively. The electron velocity is normalized on the 

electron thermal velocity Vth,e = ^D^pe = \/T e /m e and the electric field is measured in units m e vth,eW pe /e- Then, 
the dimensionless equations read: 

d t f e +vd x f e -{E + E ext )d v f e = (103) 
for the electron distribution function and 

cW = / fedv - 1 (104) 

for the electrostatic potential, <j> with E = —d x <j>- Here E cxt is an external driver added to the Vlasov equation that 
can be switched on or off during the run. The electron distribution function is discretized in space for < x < L x , 
with L x = 500 Xu the total box length, with a resolution of dx = O.IA^. The electron velocity grid ranges over 
—80 i>th,e < « < +580«th,c, with a resolution of Av = 0.0533 uth.c- Finally, periodic boundary conditions are used in 
the spatial direction. 

The plasma is initially homogeneous with waterbag electron distribution, which is modelled by the super-Gaussian 
function f e (x,v) = exp(— w 8 )/[2r(9/8)] with T(x) the Euler gamma function (34|. 

Added to the Vlasov equation Eq. ll03p external driver E ext is given by E ext (x, t) = if t < t\ or t > t%, E oxt (x, t) = 
-2A(x g /L) exp (-xQ [1 - exp (-2(i - ti))] for h < t < t 2 . Here x g = (x - x - v ph t)/L with L = 0.0625, v ph = 10, 
A = 150, h = 1 and t 2 = 3. 

The results of the Vlasov simulations of nonlinear wake wave breaking in thermal plasmas are presented in FigJTTJ 
where we show the electron phase space and electron density profile for t = 5, 6, 7. The electron momentum is 
normalized on m e u t h. and density on the ion density n . As we see in Fig. QT]a, a t time t = 5, when the electron 
velocity reaches w p h, the wake wave starts to break with the singularity corresponding to above discussed the "A- type 
breaking", which results in the narrow density spike shown in in Fig. [TT]d. The electron multistrcam region is formed 
at t > 5 as seen in Fig. [TT] b. Due to the momentum conservation the wake wave experiences a recoil leading to a 
slowing down of its propagation velocity and to a backward acceleration of the electrons in the region localized ahead 
of the wavebreaking point and to piling up the electron density, which make the electron density spike to be more 
narrow with high electron density inside (Fig. [11] b, c )- At t = 6 and 7 the electron phase space evolves into the the 
structure, which can be called "the AT- type breaking" (Fig. [TT]c, f). Later the multistream motion region becomes 
wide and the electron density maximum becomes broader. 



C. Simple model 

Consideration of Fig. [TT] showing the singularity structures formed during and after the wave breaking leads to 
the formulation of a simple model within whose framework we can explain analytically the main features seen in the 
electron density distribution. As we may sec from Fig. [TTJb the "A-typc breaking" in the phase plane, p, A, can be 
locally approximated by a superposition of two finite width stripes of parabolic and cubic form as is illustrated in Fig. 
1121 a. In other words, the waterbag distribution function is constant within the regions marked by the curves p±(X) 
given by equations 

?4/2 = A±AA/2 (105) 

in the part corresponding to the parabolic behaviour and 

P± - rp± = X - X c ± AX/2 (106) 

for the cubic part. The parameters r and X c provide the overlapping of these two stripes at large p, with A A being 
the distribution width at p = 0. 

In order to parametrize these dependences wc consider the electron motion in the frame of reference, where the 
singularity region is at rest. The parabolic stripe here can be described using an approximation of the integral of 
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FIG. 11. Results of lD-Vlasov simulations of nonlinear wake wave breaking in thermal plasmas: electron phase plane in the 
plasma described by the waterbag distribution (a,b,c) and electron density profile (d,e,f) for t — 5, 6, 7, respectively. 



motion, m e cj = m e ayo + eEX, in the vicinity of the reflection point, where p — > 0, i.e. p 2 jl = m e c(7o — 1) + eEX. 
We find that AX in Eq. (|105[) is proportional to the width Apo of the initial momentum distribution and inversely 
proportional to the electric field reflecting back the electrons in the wave breaking region: AX = (m e c/eE)(-f + ,o — 
7— ,o) ~ Apo/eE. In the laboratory frame of reference the distribution width is approximately 27 p h times narrower. 




FIG. 12. Simple model of the wave breaking: a) electron phase plane, b) electron density. 

The electron density can be calculated as the area within the p±(X) curves. Its part corresponding to the parabolic 
curves is equal to 



n e (X) = n 



2\J2eEm e c 
Ap~o 



6 X 



AX\ 



X 



AX 



e x 



AX\ 



X 



AX 



(107) 



where 6(x) is the Heaviside unit step function. The density reaches its maximum at X — AX/2 with n ei max = 
noy/m e c/ Apo- In the limit X 3> AX the electron density is inversely proportional to the square root of X, 
n e (X) ~ 1/y/X as in the case corresponding to Eq. (|102|) . In the laboratory frame of reference we have n e max = 
no ^/m e c7p h / Apo • The contribution to the electron density from the cubic part of the distribution function is pro- 
portional to the surface of the area bounded by the curves p±(X) which are the roots of equation (|106|) given by the 
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expressions 



P±(X) 



2 3 



l/3 r + 2 l/3^ 81X 2 



12r 3 - 9X± 



2/3 



6 2/3 



81Xl - 12r 3 - 9X± 



1/3 



(108) 



with X± = X — X c ± AX/2, where X is normalized on mc/eE and p measured in units of m e c. At X > Ap the 
electron density is proportional to X -2 / 3 , as in the case corresponding to Eq. (|59p. 

We see an apparent similarity between the density distribution obtained with the computer simulations, which is 
shown in Fig. [TT]f, and the density distribution given by the simple model (Fig. [12] b). 

When r = the cubic part of the distribution function develops a new breaking point and for r < it is no longer 
a single valued functions of X. At r = the contribution of the cubic part results in the electron density described by 



n e {X) 



n 
Ap 



' (x+) (x+) 1/3 + e (-x_) 1/3 - e (a_) 1/3 - {-x+) (-x+) 1/3 



(109) 



In the limit Ap — > the electron density profile for X — > is given by n e (X) ~ X 2 / 3 in accordance with the theory 
of the wave breaking in a cold plasma (see Eq. (|59p). 



D. Energy scaling of laser accelerated electrons 



Here we consider the LWFA acceleration in the above wave breaking regime when the wake field amplitude is not 
limited by the value Eap ([551) and is related via Eq. (|95p for -f m = 1 + a 2 /2 to the laser pulse amplitude as fi mas = a. 
The electron injected into the wakcficld acceleration phase can acquire the energy Q 

where the wakefield electrostatic potential is equal to tp w ~ 27moe 2 r 2 tJ with r w being the wakewave transverse size, 
which is of the order of the laser pulse waist equal to r w rj (c/uj po )\/2a. The amplitude of the laser pulse is given 
by Eq. (|85p. Using the relationship between the laser power and the amplitude (|55|) and between the wake wave 
phase velocity and the plasma density, which can be written as 7 p h = ^n CI a/no (see Ref. [3(|), we obtain for the 
accelerated electron energy 



V 



2/3 / \ 1/3 



— 1/3 

As we see, for given laser power the fast electron energy is proportional to n , i.e. the lower plasma density, 
the higher the electron energy. The electron density cannot be lower than the density determining the relativistic 
self-focusing threshold (here we do not consider the laser wakcficld excited inside a plasma waveguide, i.e. inside a 
plasma filled capillary), at which no.min = n cv V c /'P and a sa 1, i.e. the wake plasma wave is in the weakly nonlinear 
regime as required for the laser based electron-positron collider (35[, i.e. for a > 1. As the result, we obtain the 
electron energy scaling under the optimal conditions 

£ » m e c 2 f|j , (112) 

which for V = 50TW yields £ = 3GeV, and for V = 100PW gives £ = 6TcV. 

The acceleration length according to Eq. (|110p . l acc — 2r w j 2 h , in the optimal regime is given by 

A / V x 3/2 



In the case of V — 50TW one-micron wavelength laser, we have Z acc sa 5cm. 

We recall that in the limit of large laser amplitudes the energy scaling of the accelerated electrons in Eq. (|112p has 
a different dependence on the laser plasma parameters as discussed in Ref. Q and references quoted therein. 
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VII. DISCUSSIONS AND CONCLUSIONS 

In the present paper, by extending an approach formulated in Ref. fioj to the relativistic limit, we investigated the 
wave breaking of rclativistically strong Langmuir wave in thermal plasmas. As is well known, the wavebreak concept 
is meaningful only for systems which allow the hydrodynamics description because in kinetic systems with broad 
distribution functions there are always processes similar to wave breaking, such as the Landau damping in linear and 
nonlinear regimes. 

In the study of high power laser matter interaction wavebreak-like processes attract great attention in regimes where 
the wave amplitude is much larger than the distribution thermal spread in the momentum space, the most relevant 
questions being the maximal electric field, on the structure of the formed singularity and on the number of electrons 
involved in the wavebreaking. 

Using the relativistic waterbag model we showed the typical structures of singularities occurring during the wave 
breaking, we found the dependence of maximum electric field on the wave parameters, and discussed the behaviour of 
nonlinear wave in collisionless plasmas. The approach based on the warm plasma fluid model fl4| leads to the same 
scalings for the profile of the breaking waves. 

We found that in the above breaking limit the electron distribution in the nonlinear wave takes a skewed form. 
Note the somewhat similar feature in breaking water waves, when a symmetric Stokes profile [22| evolves to a skewed 
wave (see Ref. (HI). 
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